Nonlinear Structure Formation with the Environmentally Dependent Dilaton 
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We have studied the nonlinear structure formation of the environmentally dependent dilaton 
model using V-body simulations. We find that the mechanism of suppressing the scalar fifth force 
in high-density regions works very well. Within the parameter space allowed by the solar system 
tests, the dilaton model predicts small deviations of the matter power spectrum and the mass 
function from their ACDM counterparts. The importance of taking full account of the nonlinearity 
of the model is also emphasized. 



I. INTRODUCTION 

Modifying gravity on large scales is one of the plausi- 
ble ways of explaining the recent acceleration of the ex- 
pansion of the universe. So far, the construction of valid 
models of modified gravity has been fraught with difficul- 
ties. The most serious one is already present in the orig- 
inal Pauli-Fierz formulation of massive gravity [1] and 
involves the existence of a ghost in curved backgrounds 1 . 
This phenomenon seems to be generic as suggested by Os- 
trograski's theorem[2] which states that higher derivative 
theories have a Hamiltonian which is unbounded from 
below. Higher dimensional versions of modified gravity 
such as the DGP model [4] also suffer from the presence 
of a ghost in their spectrum at low energy. This problem 
is nicely avoided in f(R) models [5] which turn out to 
be equivalent to a particular type of scalar-tensor the- 
ories [6]. In these models, the compatibility with solar 
system and laboratory tests of gravity is not straightfor- 
ward and can only be achieved thanks to the so-called 
chameleon mechanism [7-12]. Indeed, the existence of a 
nearly massless scalar field on cosmological scales could 
jeopardize gravity locally. This issue is common to all 
known models of dark energy coupled to matter [13]. 
In a large class of dark energy models involving a lin- 
ear coupling to matter and a non-linear potential, the 
chameleon mechanism, whereby the scalar field mass be- 
comes dependent on the ambient environment, would be 
sufficient to hide away the dark energy field locally. Sim- 
ilarly, in models of gravity such as the DGP or Galileon 
[14, 15] theories for which a shift symmetry only allow 
for non-linearities in the scalar field kinetic terms, the 
Vainshtein mechanism [16] can be at play and prevent 
the existence of a fifth force locally. In this paper, we will 



focus on a different type of models involving a scalar field. 
These models are inspired from the string dilaton in the 
strong coupling regime [17-19]. Their gravitational va- 
lidity relies on an environmentally dependent form of the 
Damour-Polyakov mechanism [20] whereby the coupling 
to matter is driven to vanish cosmologically. Here, the 
coupling to matter is negligible in dense regions and in 
the vicinity of dense bodies. This prevents the existence 
of a fifth force in galaxies. Constraints on the parameter 
space of these models springing from local tests of gravity 
have already been obtained in [21]. Here we will study the 
cosmology of these models in the non-linear regime when 
structures form (see [22] for an analysis of the bispec- 
trum of matter distribution in this model). This requires 
large computer simulations. As a result, we have access 
to non-linear properties of the dilaton models such as the 
non-linear power spectrum or the number of dark matter 
halos for a given mass. Moreover we will be able to probe 
how much the local tests of gravity constrain large scale 
structure formation and deviations from general relativ- 
ity. We find that the dilaton models differ from GR at 
most at the level of a few percent once the local (i.e. 
solar-system) constraints have been imposed. Although 
possibly detectable in principle, observing such small de- 
viations will be challenging in the near future. 

The arrangement of this paper is as follows: in Sect. II 
we briefly review the dilaton model under study and de- 
rive the field equations in the Newtonian limit, which 
are relevant for the study of structure formation on sub- 
horizon scales and at late times. In Sect. Ill we describe 
the algorithm and code of our iV-body simulations, and 
perform relevant tests of the code. More technical details 
are given in the appendices. The numerical results and 
their analysis are summarised in Sect. IV and finally we 
conclude in Sect. V. The metric convention is (—,+,+, +) 
and wc use c = 1 unless stated otherwise. 
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II. THE ENVIRONMENTALLY DEPENDENT 
DILATON 

In this section we very briefly summarise the essen- 
tial ingredients of the environmentally dependent dilaton 
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model, which will be used for the simulations and discus- 
sions below. For more details about the model the reader 
is referred to [21]. 



A. The Model 

The dilaton model is fully specified by the following 
Einstcin-Hilbcrt action in the Einstein frame: 
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in which g is the determinant of the metric g ab , k 2 = 8irG 
with G the gravitational constant, <p is the dilaton field, 
and V(tp) its potential, which is derived from string the- 
ory in the strong coupling limit. In the matter action <S m , 
collectively represents the matter fields and A 2 (ip)g ab 
is the metric governing the geodesies of matter particles. 
In the Einstein frame, the particles feel an extra, or fifth, 
force whose strength is determined by the coupling func- 
tion /3(<p) = [\nA(ip)] where a comma denotes partial 
differentiation. The function k(tp) is given by 
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where A is a constant. Throughout this paper Latin 
indices a, b, c, ... run over 0,1,2,3 and Greek indices 
a, /3, . . . run over 1,2,3. 

Varying the action with respect to the metric g ab , we 
obtain the total energy-momentum tensor of the model, 

K 2 T ab = K, 2 A(ip)T»l - K 2 g ab V(<p) 

+k 2 {ip) [2V a ipV b ip - g ab V c ipV c (p] (3) 

where T^" is the energy-momentum tensor for fluid mat- 
ter, i.e., baryons, radiation and cold dark matter (CDM). 
Note that there is a factor A(<p) in front of T^. T£ 6 is 
not, in general, conserved but instead: 
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For pressureless dust, where = p m u a u b , u a u a = 
— 1, Eq. (4) implies that the usual continuity equation 
holds, V a (pmU a ) = 0, and hence p m is conserved. In a 
Robertson-Walker spacctime, this means that the usual 
conservation equation for matter still holds: 
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in which ' = d/dt, subscript m denotes matter, H = a/ a 
is the background expansion rate with a the scale fac- 
tor, and an overbar stands for the background value of 
a physical quantity. The gravitational field equation, or 
Einstein's equation, is given as usual: 



Gab = Rab — Ti9abR — K 2 T ab , 



(6) 



where G ab , R a b and R are respectively the Einstein ten- 
sor, Ricci tensor and Ricci scalar. 

Varying the action with respect to the scalar field ip, 
we obtain its equation of motion: 
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where T m is the trace of T"^. The energy- momentum ten- 
sor of an individual particle with mass too at position ro 
is given by 
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where r is the general spatial coordinate. Using the 
Bianchi identity we get 



-P(<p)V a <p-0(<p)<pf%, 



(9) 
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in which T^ c is the Levi-Civita connection. Clearly, if (3 = 
then this reduces to the geodesic equation in general 
relativity, as expected. 

Eqs. (6, 3, 7, 9) contain all the physics for the analysis 
below, though to implement them in TV-body simulations 
we still have to simplify them using appropriate approxi- 
mations. These will be carried out below. 

In this paper, we focus on the particular model of [21], 
which is motivated from string theory, specified by 
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where A 2 3> 1 is a parameter and ipo is the current back- 
ground value of ip. Vq is another parameter of mass di- 
mension 4. Because the potential is exponential, we are 
free to shift the value of ip so that ipo = 0. Clearly 
V c = Vbe~ Vo = Vq should be chosen carefully so that it 
can play the role of dark energy today. Both the param- 
eters A 2 and A are crucially constrained by local tests. 
In the numerical simulations, we will choose values of the 
parameters which are on the verge of the allowed parame- 
ter space in order to enhance the possible effects on large 
scales. 



B. Nonrelativistic Limit 

Eqs. (6, 3, 7, 9) are general relativistic equations. 
To implement them into iV-body simulations for large 
scale structure formation, it suffices to work in the non- 
relativistic limits, since the simulations only probe weak- 
gravity regime and small volumes compared with the cos- 
mos. 

We write the perturbed metric in the conformal New- 
tonian gauge as 
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where r, are respectively the conformal time and co- 
moving coordinate, is the metric of a 3-D Euclidean 
space, and 4>, ip respectively the Newtonian potential and 
the perturbation to the spatial curvature. For complete- 
ness, we list the expressions of the components of G a b 
in terms of the metric variables using our convention in 
Appendix A. 

Let us first look at the scalar field equation of motion 
Eq. (7). For this, we define £ such that V Q £ = fc(<£>)V a v?, 
and write V a [k(ip)V a (p\ to first order in the metric per- 
turbation variables as 

a 2 V a [*(¥>)V«d ~ -(1 - 20)r + V x £ 



2-(l- 2<£) -(<£' + 3^0 
a 



where ' = d/dr, and V x is the derivative with respect to 
the comoving coordinate x. Substituting this expression 
into Eq. (7), and removing the background equation of 
motion 
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we obtain the perturbation part of this equation: 
V x • [k(tp)V x tp] 

[%k+4%)]-%)} 

{p{(p) [A(cp)p m + iV(<p)} - V{<p)} . (16) 
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Note that in the above derivation we have dropped terms 
such as (f)', ip' and ^-<p, since we are working in the quasi- 
static limit in which the time derivative of a quantity is 
much smaller than its spatial gradient, i.e., |V X <^| ^> \4>'\- 
Using the expressions given in Appendix A, we can 
write the 00-component of the Ricci scalar as 



a 2 R° » -V x < 
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again up to first order in the perturbed metric variables. 
Similarly, 
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where T is the trace of the total energy-momentum ten- 
sor. Then the 00-component of the Einstein equation 
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with the background part, i.e., the Raychaudhuri equa- 
tion, 
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removed, can be written as 
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- A{ip)p B 
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where we have defined $ = a<j) for convenience. 

Finally, for the equation of motion of matter particles, 
Eq. (9), using the relationship between physical coordi- 
nates r and comoving distance x, we can rewrite it as 
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Defining the conjugate momentum to x as p = a 2 x, this 
equation could be decomposed as 



dx 
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dp 
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Note that there are two components of the fifth force, as 
discussed in [25] 

Eqs. (16, 19, 21, 22) are all that we need to put into the 
iV-body simulation code to study structure formation in 
the nonlinear regime. Before that we have to discretise 
these equations and write them using code units, so that 
they can be applied on a mesh with finite grid size. These 
lengthy expressions are given in Appendix B, where we 
also discuss the subtleties in the numerical implementa- 
tion. 



III. THE iV-BODY SIMULATIONS 

In this section we briefly describe the algorithm and 
model specifications of the TV-body simulations we have 
performed. We also give results for the tests of the code, 
which show that the scalar-field solver works quite well. 

A. Outline of the Simulation Algorithm 

For our simulations we have used a modified version of 
the publicly-available A-body code MLAPM [27]. The mod- 
ifications we have made follow the detailed prescription 
of Rcf. [25] , and here we only give a brief description. 

The MLAPM code has two sets of meshes: the first in- 
cludes a series of increasingly refined regular meshes cov- 
ering the whole cubic simulation box, with respectively 
4, 8, 16, ■ ■ • , N<i cells on each side, where A^ is the size of 
the domain grid, which is the most refined of these reg- 
ular meshes. This set of meshes arc needed to solve the 
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Poisson equation using multigrid method or fast Fourier 
transform (for the latter only the domain grid is neces- 
sary). When the particle density in a cell exceeds a pre- 
defined threshold, the cell is further refined into eight 
equally sized cubic cells; the refinement is done on a 
cell-by-cell basis and the resulting refinement could have 
arbitrary shape which matches the true equal-density 
contours of the matter distribution. This second set of 
meshes are used to solve the Poisson equation using the 
linear Gauss-Seidel relaxation scheme. 

The dilaton field is the most important ingredient in 
the model studied here, and we have to solve it to obtain 
detailed information about the fifth force. In our TV-body 
code, we have added a new scalar field solver which is 
based on Eqs. (B13, B14, B15, B16). It uses a nonlin- 
ear Gauss-Seidel scheme for the relaxation iteration and 
the same criterion for convergence as the default Pois- 
son solver in MLAPM. But it uses V-cycle [28] instead of 
the self-adaptive scheme in arranging the Gauss-Seidel 
iterations. 

The value of u (see definition in Appendix B) solved in 
this way is then used to calculate the total energy den- 
sity including that of the scalar field, and this completes 
the computation of the source term to the modified Pois- 
son equation. The latter is then solved using fast Fourier 
transform on the domain grid and Gauss-Seidel relax- 
ation on refinements, according to Eq. (B18). 

With the gravitational potential $ and the scalar field 
u at hand, we can use Eq. (B'20) to evaluate the total force 
on the particles and update their momenta/velocities. 
Then Eq. (B19) is used to advance the particles in space. 

For more details about the implementation see [25]. 



B. Simulation Details 

The physical parameters we use in the simulations are 
as follows: the present dark-energy fractional energy den- 
sity tt A = 0.743 and Q m = 0.257, H = 71.9 km/s/Mpc, 
n s = 0.963 and as = 0.769. We use two sets of simula- 
tion box which have sizes of 32ft, -1 Mpc and 64ft -1 Mpc 
respectively, in which ft = Hq/(100 km/s/Mpc). We sim- 
ulate four models, with parameters (A 2} A) = (4x 10 6 , 2), 
(4 x 10 5 , 10), (2 x 10 5 , 100) and (2 x 10 6 , 30). These pa- 
rameters are chosen so that they predict local fifth forces 
which are allowed by current experiments and observa- 
tions 2 . In all those simulations, the particle number is 
256 3 , so that the mass resolution is 1.114 x 10 9 ft _1 Mq 
for the 64ft,- 1 Mpc simulations and 1.393 x 10 8 ft _1 Mq 
for the 32ft -1 Mpc simulations. The domain grid is a 
128 x 128 x 128 cubic and the finest refined grids have 




2 The values are taken from near the boundary of the allowed 
region in the parameter space in Fig. 1 of [21], As a result we 
expect that they should give us the biggest effect on large-scale 
structure while satisfying constraint from local experiments. 



x (h'' Mpc) 

FIG. 1. A first test of the scalar field solver. For this test we 
use a simulation box of 256ft _1 Mpc on each side, and set the 
density field to be homogeneous in the box. The exact value of 
13, /3 t h, is known analytically. The differences between /3 t h the 
initial guess of /3 in the grid cells along the s-axis are shown as 
symbols, while that between /3 t h and the f3 after relaxation are 
shown as the continuous curve. Clearly the relaxation works 
accurately. 



16384 cells on each side, corresponding to a force resolu- 
tion of about 12ft -1 kpc and 6ft -1 kpc respectively for 
the two sets of simulations. The force resolution deter- 
mines the smallest scale on which the numerical results 
are reliable. We have also run a ACDM simulation with 
the same physical parameters. 

Our simulations are purely iV-body, which means that 
baryonic physics has not been included in the numerical 
code. We use the same initial conditions for the dilaton 
and the ACDM simulations, because before the initial 
redshift Zi = 49 the fifth force is strongly suppressed 
so that the effect of the dilaton on the matter power 
spectrum is negligible. 



C. Code Tests 

Before displaying the numerical results from the N- 
body simulations, we show some evidence that our nu- 
merical procedure works correctly. As our modification 
to the default MLAPM code is only in the scalar field part, 
we focus on tests of the scalar field solver and the fifth 
force only. 

The scalar field solver uses the nonlinear Gauss-Seidel 
relaxation scheme to compute (3, and an indicator that 
it works is to show that, given the initial guess of the 
solution that is very different from the true solution, the 
relaxation could produce the latter within a reasonable 
number of iterations. Consider a simulation box with ho- 
mogeneous density, then the true solution to /3, /3 tn , could 
be calculated analytically. We therefore make an initial 
guess for j3 which is randomly scattered around /3th and 
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1 10 
x (h"' Mpc) 

FIG. 2. A second test of the scalar field solver. For this test 
we use a simulation box of 128/i _1 Mpc on each side, and set 
the density field to be the equivalent of having a point mass at 
x — y — z — and zero otherwise. Far from the point mass, 
the solution to ft can be approximately solved analytically 
(the continuous curve). The symbols show the results for ft 
from the numerical code. The two show good agreement in a 
wide range of x. 



let the scalar solver try to recover /3 t h- In Fig- 1 we have 
shown |/3 — /3 t h | before (symbols) and after (curve) the 
relaxation: as can be seen there, before the relaxation 
the difference between the initial guess (3 and /3th is of 
order 0.01, while after the relaxation it reduces to 10~ 6 . 
Note that 10~ 6 corresponds to the error caused by using 
floating-point numbers, and as a result this shows that 
the scalar solver works accurately. 

As a second test of the scalar field solver, consider hav- 
ing a point mass at the origin x = y = z = and the vac- 
uum density otherwise. This could be achieved by filling 
the densities in the cells of the simulation grid according 
to [29] 

Pc = 10- 4 /V d 3 (23) 



Using the code units (see Appendix B), this can be writ- 
ten as 

V 2 <5/3 » m 2 ca S(3 (26) 

with 

2 (M ) 2 on 3 3^(1 - 2P) + 2A- 2 

mta = 5 — 3i2Aa°^42 5 — op. (27) 

«c 2 ( 3/3 2 + A -2) 2 

Here B is the box size of the simulation box, and c is 
the speed of light, which we have restored to make the 
dimension explicit. The analytic solution is thus 

Sj3( r ) = - e ~ m " ltr (28) 
r 

where C is some constant determined by the value of the 
point mass and r the distance from the origin. Because 
C is unknown, we fix its value by requiring that Eq. (28) 
be equal to the numerical solution at r = 10ft -1 Mpc. 

The normalised analytical solution to /3 is shown as the 
continuous curve in Fig. 2, while the numerical solutions 
are shown as symbols. We see that the two agree over a 
wide range of r (note that \S/3\ changes by several orders 
of magnitude) . Note that when r is small the agreement is 
not perfect, because linearisation does not work very well 
near the high density region; meanwhile, for very big r the 
value of \Sj3\ drops below O(10 -6 ) and numerical error 
due to using floating-point numbers becomes important. 

In summary, Figs. 1 and 2 show that our scalar solver 
works well. Below, we also show that the fifth force agrees 
with analytic approximations in certain regimes. 

IV. NUMERICAL RESULTS 

In this section, we shall present our simulation results, 
including the snapshots, the matter power spectrum and 
the halo mass function. 



A. Snapshots 



for the cell with i = j = k = 0, in which Nd is the number 
of cells on each side of the domain grid, and p c = 10 -4 
for all other cells. 

Outside the particle it is the vacuum, in which the 
scalar field equation of motion can be approximately lin- 
earised as 



V 2 J(3 w 8irGV A 2 a 2 3/? ^ 



2/3) + 2\- 



(W 2 



X- 



-6/3, (24) 



where we remind the reader that /3 = A 2 <p (see eq. 10 
with ipo = 0), and 5/3 = (3 — /3. (3 is the background value 
of /3, which can be analytically calculated as [21] 



(25) 



a. 



40 A a 3 



As we have seen, in the dilaton model /3 and thus the 
fifth force is suppressed in high density regions. In this 
subsection we demonstrate these qualitative features us- 
ing some snapshots. 

Fig. 3 shows the comparison of the magnitudes of the 
fifth force and gravity for the four models we have con- 
sidered, at three different output times a = 0.2,0.5 and 
1.0. For this, we pick out a thin slice from the middle 
of the z = I6/1 -1 Mpc simulation box, and compute the 
fifth force and gravity on the particles within that slice. 

At early times, the density is high everywhere and we 
expect the fifth force on all particles in all the four models 
to be strongly suppressed, and this is confirmed by the 
first row of Fig. 3, which shows that the fifth force is much 
weaker than gravity. Note that the degree of suppression 
of the fifth force is dependent on the value of A 2 : the 
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FIG. 3. The magnitude of the fifth force (the vertical axis) versus that of gravity (the horizontal axis), for the particles (black 
points) selected from a thin slice of the 32/i _1 Mpc simulation box. We show this for the four models we have simulated and at 
three different output times, as given in the subtitle of each panel. Note that both forces are expressed using the internal unit 
(see Appendix B), which is Hq/B times the physical force unit. 



larger A 2 is, the more the fifth force is suppressed. Also, 
the fifth force is weaker in higher density regions (where 
gravity is stronger) than in lower density regions (where 
gravity is weaker), showing a strong dependence on the 
environment. 

As the Universe expands, the overall density decreases 
and the fifth force becomes stronger, which could be seen 
in the lower rows of Fig. 3. If there is no suppression on 
the fifth force, then its strength should be 

a = WTx^ (29) 

times that of gravity [21]. For comparison, in the low- 
est row (a = 1) we have over-plotted a times gravity as 
continuous curves. We can see that in the model with 
smaller A 2 (the middle two columns) Eq. (29) gives a 
relatively good description of the fifth force at least in 
some regions. But for the models with big A 2 (the first 
and fourth columns) the fifth force is strongly suppressed 
even today. 

Since (3 determines the strength of the fifth force 
[cf. Eq. 9], we are also interested in it. Fig. 4 shows the 



values of /? as a function of position in the same slices as 
Fig. 3. As expected, at very early times (o = 0.2) j3 <C 1 
because the fifth force is strongly suppressed. As the Uni- 
verse expands, (3 increases (the colour on the points be- 
comes blue rather than black), but in the high density 
regions /3 remains very small. Also, for the models with 
big A 2 (the first and fourth columns) the values of (3 in 
high and low density regions tend to have stronger con- 
trast, showing stronger environment-dependence. This is 
clearer in the third row, which shows the result at a = 1.0. 
This is consistent with what we have seen in Fig. 3. 



B. Matter Power Spectrum 

The nonlinear matter power spectrum is an important 
structure formation observable and could be used to dis- 
tinguish amongst different structure formation scenarios. 
In Ref. [21] it has been shown that the growth rate of 
linear matter density perturbations in the dilaton model 
differs from that of ACDM only slightly, and therefore the 
dilaton model (with its parameters constrained by solar 
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FIG. 4. (Colour Online) The colour scale plot of the value of /3 as a function of coordinates x, y, for the same thin slice of 
the 32h~ 1 Mpc simulation box as in Fig. 3. We show this for the four models we have simulated and at three different output 
times, as given in the subtitle of each panel. Each point represents a particle, and the colour of the point depends on the value 
of jS at the position of that particle: for all the panels the lightest colour (white) denotes p — 0.3 and the darkest colour (black) 
denotes j3 = 10 -7 ; the blue colour is interpolated linearly between these two extremes 



system tests) does not deviate at more than the percent 
level from ACDM in practice. On the other hand, the fifth 
force in the dilaton model has a finite range and is ex- 
pected to only take effect on the scales of galaxy clusters 
(~ O(Mpc)) and smaller, which already fall into the non- 
linear regime. We are therefore interested in seeing how 
the fifth force affects the growth of density perturbations 
on these scales. 

Fig. 5 displays the fractional difference of the dilaton 
nonlinear matter power spectrum from that of the ACDM 
model, defined as (-P(fc) - -PacdmW) / -Pacdm From 
this we can see that the difference is strongly suppressed 
even on small scales where the fifth force is expected to 
take effect. This is different from the linear perturbation 
prediction of [21] (c.f. Fig. 3 therein), which shows that 
the growth rate of density perturbation on small scales is 
significantly higher than that on large scales. The reason 
for this is that by linearising the scalar field equation, 
the nonlincarity of the dilaton model, which is the very 
mechanism that suppresses the fifth force in high density 
regions, is artificially removed (at least partially), and the 



strength of the fifth force is determined by the average, 
instead of the local, matter density. In contrast, the N- 
body simulation overcomes this problem by taking full 
account of the suppression of the fifth force. 

The results indicate that it is even more difficult to 
use the nonlinear matter power spectrum to constrain 
the dilaton model or distinguish it from ACDM as the 
differences to the ACDM power spectrum are only a few 
percent on very small length scales at late times. 



C. Mass Function 

The halo mass function is another key structure forma- 
tion observable. It is defined to be the number density of 
dark matter halos within a given mass range. Clearly, in 
case of a fifth force which could boost the clustering of 
matter, we expect more halos to form. In Fig. 6 we have 
shown the mass functions of the dilaton models compared 
with that of ACDM, at z — 0. Although the dilaton mod- 
els do have higher mass functions than ACDM, especially 
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FIG. 5. (Colour Online) The fractional difference between the dilaton and ACDM nonlinear matter power spectra, which is 
defined to be (P(k) — Pacdm(&)) / -Pacdm (fc) ■ The black, green, pink and purple curves are respectively for the models with 
(A 2 ,X) = (4 x 10 6 ,2), (4 x 10 5 ,10), (2 x 10 5 ,100) and (2 x 10 6 ,30). The four panels (upper-left, upper-right, lower-left and 
lower-right) are results at a = 0.3, 0.5, 0.7 and 1.0. 



for small halos which generally live in low-density regions 
where the fifth force is less suppressed, the differences are 
again very small, making all these models hard to distin- 
guish in practice at present and a challenge for future 
surveys. 



D. Halo Profile of /3 

In Fig. 7 we show the profile of (3 inside the dark mat- 
ter halos, which are assumed to be spherical. Because /3 
characterises the strength of the fifth force, this can also 
provide information about the fifth force in halos. We 
have selected three halos with different masses (respec- 
tively 3.5 x 10 14 , 6.0 x 10 13 and 1.6 x 10 13 solar mass) to 
check the results. 

As can be seen from Fig. 7, j3 (and thus the strength 
of the fifth force) increases from inner to outer regions of 
the halos as the matter density is highest and the fifth 
force most severely suppressed in the central region. Fur- 



thermore, the fifth force is stronger for smaller halos, be- 
cause those generally reside in low-density regions where 
the fifth force is less suppressed. However, for all these 
selected halos (3 is at most ~ 0(1 CP 2 ) and typically less 
than ~ O(10~ 3 ) except near the halo edge, which mean 
that well inside the halos (such as where the solar system 
is) the fifth force is much weaker than gravity and has 
negligible effects (and we have not even included baryons 
in the simulations, which arc generally much denser than 
dark matter in galaxies) . 



The results are consistent with what we have seen in 
the nonlinear matter power spectra and mass functions, 
all of them showing that the fifth force has little influence 
in the structure formation of the dilaton models (as long 
as solar system tests are passed). 
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FIG. 6. (Colour Online) The mass functions for the ACDM 
model (black solid curve), and the four dilation models with 
(A 2 , A) = (4x 10 6 ,2), (4 x 10 s , 10), (2 x 10 5 ,100) and (2 x 
10 6 , 30) (the coloured curves). All the curves are very close. 



which reside in voids, where the overall density is low and 
the fifth force could be as strong as gravity. However, the 
spatial and mass resolutions of our simulations do not 
allow a detailed analysis of this. 

Note that in the simulations of this work we have only 
included dark matter but not baryons. However, as long 
as the scalar field has a uniform coupling to different mat- 
ter species, we expect that all the results will qualitatively 
remain. In particular, in the inner regions of the halos, 
baryon density is much higher than that of dark matter, 
which could further suppress the fifth force compared to 
what we have seen in our simulations. 

In conclusion, while obscrvationally hard to distinguish 
from the ACDM model, the environmentally dependent 
dilaton model is a very effective way to shield the dilaton 
from observations. 



SUMMARY AND CONCLUSION 
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Appendix A: Useful Expressions 



Up to first order in the perturbed metric variables (j), ip, 
the nonzero components of the symmetric Levi-Civita 
connection arc 
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FIG. 7. The profiles for /3(r) in some chosen halos. The diamond, triangle and box represent results for three halos with masses 
equal to 3.5 x 10 14 , 6.0 x 10 13 and 1.6 x 10 13 h Mq respectively. The horizontal axis is the distance from the halo centre, in 
unit of h~ kpc. 



The components of the Ricci tensor and Ricci scalar up 
to first order in <f>, ip are then easy to compute as 
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Appendix B: Discretisation of Equations 

To implement the nonrclativistic equations into our 
numerical code, we have to rewrite them using code units, 
which are given by 



x 

x c Pc 



-, t c — tH t 



0, 



<I>r 



H B 

^-a, Pc = ^— , V c = BV X (Bl) 



in which a subscript c denotes code unit, .B is the size of 
the simulation box and Hq = lOO/i km/s/Mpc. In what 
follows we shall write V = V c for simplicity. 



1. Scalar Field Equation of Motion 

Recall that in our model we have chosen ipo = such 
that /3(<p) = A-2<p. It is then the same to solve for ip or 
to solve for /3. As A% ^> 1, j3 ^> <p and so we shall solve 
for (3 rather than ip to reduce possible numerical errors. 
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The equation of motion for (3 could be obtained simply 
by multiplying that for tp by A 2 : 

V x • [k(p)V x p] 
4irGA 2 a 2 



4nGA 2 a 2 



[P{A(P)p m +AV(P)]-V(P)} 
[p {A(p)p m + W{fi)] - V(p)} ,(B2) 



where we have used P instead of ip as the variable. 

As discussed in [21], j3 characterises the strength of the 
fifth force. In high density environments, f3 -c 1 so that 
the fifth force is too weak to be measured; in low den- 
sity regions, however, we have P ~ 0.23 today, indicating 
that the fifth force is roughly as strong as gravity. Fur- 
thermore, Eq. (B2) does not say anything about the sign 
of P. 

Obviously, because /3 ranges from O (10~ 6 ) to 0(1), 
using P directly in the numerical code could easily cause 
big numerical errors in the regions where P is small. One 
alternative is to use ln(/3) as a new variable, but this 



does not necessarily work because P might be negative. 
Therefore, in this work, we shall use a different variable 
u = with n being some odd positive integer, as 

the redefined scalar field. More explicitly, we shall adopt 
n = 9 which guarantees that u ~ 0(0.1 — 1), i.e., u spans 
a much smaller range than p, making it easier to control 
numerical errors. Furthermore, n being odd makes sure 
that u is never undefined even if P < 0. 
In terms of u, we have 



k(u) = \Jiu 2n + A- 2 , 

,2n 



A{u) = 1 + 
V{u) 
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and using the code units defined above, we could rewrite 
the scalar equation of motion Eq. (B2) as 
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where we have defined 



b(u) = nu n - V3u 2n + A- 



and P is the background value of P, which can be com- 
puted as [21] 



where a subscript j means that the quantity is evalu- 
ated on the j-th point. The generalisation to the three 

(B8) 

dimensional case is straightforward. 

The factor b in V • (6Vu) makes this a standard variable 
coefficient problem. We need also to discretise 6, and do 
it in this way (again for one dimension) [23]: 



P = 



(B9) 



The full equation for u, Eq. (B7), contains the quan- 
tity V • (6Vu). To discretise it, we shall assume that the 
discretisation is performed on a grid with grid spacing h. 
We shall require second order precision which is the same 
as the default Poisson solver in MLAPM, and then Vm in 
one dimension can be written as 



V h Ui 
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in which bj±i = \(bj +bj±i). Generalising this to three 
(BIO) dimensions, we have 
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Then the discrete version of Eq. (B7) is 

L h (u^k) = 0, (B13) 
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Then the Newton- Gauss- Seidel iteration says that we can 
obtain a new (and usually more accurate) solution of u, 
u ?Tfc> usm S our knowledge about the old (and less acu- 



rate) solution 
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The old solution will be replaced with the new one once 
the latter is ready, using a red-black Gauss-Scidcl sweep- 
ing scheme. Note that 
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where we have defined 
_ db{u) 



d(u) 



du 

3n(2n - l)u 3n ~ 2 + n(n - l)A~ 2 u 
V3w 2 " + A- 2 " 



2„,n,-2 



.(B17) 



ground value because we expect that any perturbations 
should be small then. For subsequent time steps we can 
use either the solution at the last time step or some an- 
alytical approximated solution as the initial guess. 



2. Poisson Equation 



In principle, if we start from some high redshift, then 
the initial guess of Uij.fc could be chosen as the back- 



In terms of the newly-defined scalar field u and using 
the code units, the modified Poisson equation becomes 
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The discretisation of V 2 <E> C is straightforward and will 
not be presented here. 
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3. Particle Equation of Motion 

Using the code units, Eq. (21) could be easily rewritten 

as 



dx c 

dt c 



(B19) 
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Similarly, Eq. (22) becomes 

-P = V$ c hhH -Vw 
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